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Analysis of a Reduced-Order Approximate Deconvolution 
Model and its interpretation as a Navier-Stokes-Voigt 

regularization 
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Abstract 

We study mathematical and physical properties of a family of recently introduced, reduced- 
order approximate deconvolution models. We first show a connection between these models and 
the NS-Voigt model, and that NS-Voigt can be re-derived in the approximate deconvolution 
framework. We then study the energy balance and spectra of the model, and provide results of 
some turbulent flow computations that backs up the theory. Analysis of global attractors for 
the model is also provided, as is a detailed analysis of the Voigt model’s treatment of pulsatile 
flow. 


1 Introduction 

Over the past several decades, many large eddy simulation (LES) models have been introduced 
and tested, with varying degrees of success; see e.g. mug for extensive overviews. Each of these 
models has its share of good and bad properties, and with each generation of models, the amount of 
“good” tends to increase while the amount of “bad” tends to decrease. Two attractive models of the 
current generation are approximate deconvolution models (ADMs) [U :2] , and Voigt regularization 
models [13 Ell, both of which have recently seen significant interest from both mathematicians 
and engineers. These models are known to be well-posed hi H3 na El], have attractive physical 
properties (such as energy and helicity conservation and cascades [3H [35] ), and have performed well 
in (at least some) numerical simulations !35) iL5], including when applied to the Euler equations for 
ideal fluids m- 

However, both models have their downsides: It is unknown how to efficiently and stably imple¬ 
ment the usual ADM with standard finite element methods, and the Voigt regularization model can 
lack accuracy, likely due to its low consistency error with both the Navier-Stokes equations (NSE) 

Ut + {u ■ V) u — vAu + Vp = /, (1-1) 

V • u = 0, (1.2) 

and the spatially filtered NSE E3IS2- Here u denotes the fluid velocity and p the kinematic pressure. 

In the recent papers [23 SI], a new system was proposed as an alteration of the usual ADM, but 
which is amenable to unconditionally stable and efficient implementation in a finite element context. 
Here, we add further insight in the derivation, from an abstract point of view, see Section This 
model, which we study herein, is referred to as the reduced-order ADM (RADM), and is given by 

v t — a 2 Av t + Dv ■ VDv + V<? — vADv = /, (1.3) 

V-Dv = 0, (1.4) 
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where the vector v represents an averaged velocity, q a modified pressure, while D is an approximate 
deconvolution operator meant to represent an approximate inverse to the Helmholtz filter F, with 
Fcj) := (f> satisfying the filter equation: 

— a 2 Acf) + </> = </>. 


The evolution system (1.3)-(1.4) requires an initial condition for velocity, and both the filter and 
evolution systems require appropriate boundary conditions. By restricting to the space-periodic 
setting (as usual for these problems, with only a very few remarkable exceptions), we can assume 
D to be linear, positive, self-adjoint, to commute with F, and moreover that there exists do, d± € K. 
such that 

1 < do < ||D|| < d\ < oo, 

where ||Z)|| denotes the operator norm of D. These assumptions are true for most common types of 
approximate deconvolution defined with the Helmholtz filtering operation, including van Cittert US], 
multiscale m, Tikhonov, and Tikhonov-Lavrentiev [33], see also the review in |6j. While we will 
often be general in our use of D, in some circumstances such as computations and energy spectra 
calculations, it is necessary to select a specific operator. In these cases, we will use the most common 
type of approximate deconvolution, van Cittert, and we will denote D — Dn, where N denotes the 
order of deconvolution. This operator is defined in detail in Section [2] 

The RADM system (1.3)-(1.4) was derived from the fourth-order version of the ADM by making 
the approximation D ss F~ x in the viscous term, which reduced the fourth order system to second 
order. This system was proved to be well posed in [41| and of high order consistency with the spatially 
filtered NSE. In [26], with D chosen to be van Cittert approximate deconvolution, an efficient finite 
element implementation of the RADM performed very well on the Re T = 180 turbulent channel flow 
benchmark test, as well as in a 2D benchmark flow problem with a contraction and two outlets. 

Plan of the Paper. It is the purpose of this paper to study some fundamental mathematical 
and physical properties of the model (1.3)-(1.4). After providing mathematical preliminaries in 
Section 


Section 


we show a strong connection between the RADM (1.3) and the NS-Voigt model (3.5) in 


In fact, we show that the RADM can be considered as a NS-Voigt-deconvolution model. 
That is, if the regularization operator P -1 is replaced with D~ 1 F ~ 1 , then the RADM is recovered 
after a change of variables. This is interesting and important since a significant amount of theory 
has been developed for the NS-Voigt model by Titi and coworkers in e.g. [251 S3 S2] as well as 
in [7]. We also show that NS-Voigt can be derived using the approximate deconvolution approach 
originally developed by Stolz and Adams d 1 SI S5]. In this respect we establish a strong tie 
between Voigt-type models and approximate deconvolution type models in general. 

An energy analysis of the RADM is given in Section 4. This includes an energy balance and a 
Kraichnan-type analysis of the energy spectra of the model. Some numerical results are also given 
here, to back up the theory. Our results reveal that the RADM is more computable than the NSE 
and has spectral scaling similar to the Leray-cc model [18], which matches that of true fluid flow on 
larger scales in the inertial range, but smaller scales are transferred more quickly to even smaller 
scales. Moreover, our analysis shows that increasing the deconvolution order N shortens the inertial 
range and creates less energy in each scale. 

Section 5 considers long-time behavior of the model, and proves existence of smooth global 
attractors. Section 6 considers the pseudo-parabolic nature of the system, and presents a study 
of some analytical properties and the model’s treatment of pulsatile flow, which is an interesting 
example because it allows for exact analytical solutions. Finally, conclusions are drawn in the final 
section, and prospective future work is discussed. 


2 Preliminaries 

We present here notation and mathematical preliminaries to make the analysis of later sections 
smoother. As is typical with theoretical analyses of fluid models, we restrict our attention herein 
to problems involving periodic flows on the box D = (0, 2nL) 3 . We stress that all machinery 
requires the space-periodic setting to have the commutation between the deconvolution operator 
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and general differential operators. The application of ADM methods to motion bounded by rigid 
walls is problematic and the theory is still missing. Some results in special geometries can be found 
in 0051. In this framework, NS-Voigt represents one of the few models which is naturally set in a 
bounded domain with Diriclilet conditions. 

Furthermore, we consider the system (1.3)-(1.4) to be non-dimensionalized by introducing L, 
U, and L/U for the characteristic length, velocity, and time scales. Thus, the viscosity v > 0 can 
be considered dimensionless and representing the inverse of the Reynolds number of the flow. To 
simplify the notation, in the theorems we will suppose without loss of generality that L = 1. 

We characterize the divergence-free spaces by using Fourier Series on the 3D torus. We denote 
by (ei,e 2 ,e 3 ) the orthonormal basis of M 3 , and by x := (xi,£ 2 ,£ 3 ) £ M 3 the standard point in 
M 3 . Let T be the torus defined by T := R 3 /Z 3 . We use || • || and (•,•) to denote the L 2 ( T) norm 
and inner product, and we impose the zero mean-value condition on velocity, pressure, and external 
force. We will use the customary Lebesgue L p ( T) and Sobolev W k ’ p (T) spaces, and for simplicity 
(if not explicitly needed) we do not distinguish between scalar and vector valued real functions. 

We define, for any exponent s > 0, 


V s := jw:T^R 3 , w£ W S ’ 2 (T) 3 , V ■ w = 0, J w(x)dx = 0 

where IF S,2 (T) 3 := [1F S,2 (T)] 3 and if 0 < s < 1 the condition V • w = 0 must be understood in a 
weak sense. For w £ V s , we can expand the velocity field with Fourier series as follows 

w(x) = ^ Wke lk ' x , where k = (ki, & 2 , k 3 ) £ Z 3 is the wave-number vector, 

|fc|>i 

and the Fourier coefficients are given by C 3 9 Wk ■= |yj- f T w(x)e~ lk ' x dx, with |T| the Lebesgue 
measure of T. If |fc| := \/\ki\ 2 + Ifel 2 + \k 3 \ 2 , then the V s norm is defined by 

IMIv« : = H \ k \ 2s \™k\ 2 , 

l*l>i 

where, as above, ||iu||jjo := ||iu||. The inner products associated with these norms are 

( w,v) V s = ^2 \k\ 2s w k -v k . 

W>1 

We can finally characterize V s C fF s,2 (T) 3 as follows: 

V s := jw(x) = ^ w k e lk x : ^ |fc| 2s |u; fc | 2 < 00 , k ■ w k = 0, W- k = w k 

\k\>\ |fc|>l 

Recall that due to the zero mean value, the Poincare inequality holds in V 1 : There exists A > 0 
such that for every (j> £ V 1 , 

M < A" 1 / 2 ||V0|| . 

This inequality will be used extensively in our analysis, and allows us to use the more convenient 
gradient norm in place of the V 1 norm. 


2.1 Approximate deconvolution 

We assume the deconvolution operator D is self adjoint and commutes with both A and F. These 
assumptions are known to be satisfied for van Cittert approximate deconvolution |191 i,35j , multiscale 
deconvolution |2()| . Tikhonov, and Tikhonov-Lavrentiev [33]. See also the reviews in ling for more 
results about these operators. 

We denote the filtering operator A := F -1 = (—a 2 A + J), and we introduce now a family of 
abstract deconvolution operators {Dn}n, which are of order zero, and such that 

D 0 = I and Dm —> A as N —> + 00 . 
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We note that van Cittert approximate deconvolution, 

N 

D N = J2( I ~ F ) n , 

n—0 

forms such a family, and the other types of deconvolution mentioned above can also be written to 
satisfy the above limiting property. Although we will use a general deconvolution operator D when 
possible, in some analysis (e.g. energy spectra analysis) and for computations, it is necessary to 
specify a specific operator. In these cases, we will use van Cittert approximate deconvolution, as it 
is one of the most popular in the ADM community. The basic properties we will use of the operator 
D = Djy are summarized in the following lemma, see m- 

Lemma 2.1. For each N £ N the operator Djv : V s —* V s is self-adjoint, it commutes with 
differentiation, and the following properties hold true: If Dn denotes the symbol of the operator Dn, 
then 


1 < D N (k) <N+1 


lim Dj^(k) = N + 1 

|/c|—>-+oo 


D N (k)<(l+a 2 \k\ 2 ) 


Vfc g Z 3 , 
for large |fc|, 
for fixed a > 0, 
Vfc e Z 3 , a > 0. 


This lemma shows that in the specific case of practical use, we can assume that the upper and 
lower bounds for Dn are do = 1 and di = N + 1, which reflects in similar bounds for D^. 


2.2 Known theory for RADM solutions 

We define now the notion of weak solution for the RADM. 

Definition 2.1. We say that v € L°°(0, T; IV 1 ) is a weak solution to (1.3)-(1.4) if 
d 


dt 


[(u, <j>) + a 2 (Vv, V0)] + (Dv ■ VDv, </>) + v(VDv, Vf) = (/, <j), 


for all (j) G V 1 , and also provided that Vt G L 2 (0,T;V - 1 ) and the following energy identity holds 
true: 

Vfl 1/2 r(T) \ D 1 / 2 v{T) 2 ^+"f o ||VDu|| 2 dt 

= \(o? VD 1/2 u(0) 2 + D 1 ' 2 v{ 0) 2 ^)+j o ( f,Dv)dt . (2.1) 


We have the following existence and uniqueness theorem, which is based on a standard Galerkin 
approach, together with the energy estimate (2.1l below, see HU Thm. 2.1], 

Theorem 2.1. Let f £ L 2 (0, T; ff -1 (n)) and vq £ V 1 be given. Then, there exists a unique weak 
solution to (1.3)-(1.4). 

The proof of this result is standard and follows the same guidelines of the results in dun], 
with the key point being the energy identity (2.1|. From [41], with the assumptions on the data, 
we have that the solution v £ L°°(0, T; V 1 ), and so all manipulations in this proof are justified. 
Multiplying (|1.3[) by Dv and integrating over the physical domain gives 


a 2 (Vr t , \7Dv) + (v t , Dv) + (Dv ■ VDv , Dv) — (g, V ■ Dv) + v || VDv\\ = (/, Dv), 
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and then using a vector identity on the first two terms, and that V • Dv = 0, both the pressure and 
nonlinear terms vanish, leaving us with 


1 d_ 

2 dt 



VD 1/2 v 


2 

+ 


D x ' 2 v 


v\\VDv\\ 2 = (f,Dv). 


Integrating in time over [0, T] provides (2.1). 

This equality shows that the Galerkin approximate solution v m = Yji<\k\<m x is bounded 
uniformly in L°°(0,T; V 1 ). Then, by making use of the comparison argument one can easily prove 
that 

v™ - a 2 Av™ £ L 2 (0, T; V" 1 ). 


The Lax-Milgram lemma set in the space V 1 implies that v™ £ L 2 (0,T; V 1 ). In particular, this 
proves that one can extract from the Galerkin approximate functions a converging sub-sequence, and 
moreover shows also that the solution v £ C([0, T]; V 1 ), which is a uniqueness class. The regularity 
of the time derivative allows to show that an energy equality (instead of an inequality as for the 
NSE) holds true. 


In the case of the problem without viscosity, i.e. the RADM-Euler equations 


v t — a 2 Av t + Dv ■ X7Dv + X7q = /, 
V ■ Dv = 0, 


( 2 . 2 ) 

(2.3) 


by adapting techniques from 13 na and by using the estimates on the operator D from Lemma |2.1| 
it is easy to infer the following result (the same applies also to the Euler-Voigt-deconvolution model, 
as discussed in the next section) 


Theorem 2.2. Let T > 0, m > 1, and v 0 £ V m . Let f £ C([—T, T]; W m_1,6 / 5 ) with V • / = 0. 
Then, there exists a unique solution v of the RADM Euler equations (2.2 1-( [2)3| which belongs to 
C 1 {[-T,T]-V m ). Moreover, 


sup ||v(i)||*™» < C(a, ||u 0 ||^m, sup ||/(f)|| wm -i,6/5,T,d 0 ,di)- 
te[-T,T] —T <t<T 


3 Remarks on the derivation of the model and a connection 
to NS-Voigt 

Interestingly, despite its derivation as an altered/reduced ADM, the proof for well-posedness of the 
system in m uses ideas of the well-posedness proof for NS-Voigt from [7j. This suggests a connection 
between the two models, and we show now that a NS-Voigt deconvolution model can be derived 


from (1.3)-(1.41. Writing 


v t - a 2 Av t = F~ l v t = F~ L D~ i Dv t = D~ l F~\Dv) 


and setting w = Dv yields 


D 1 F l w t + w ■ Vw + Vg — uAw = /, 

V • w = 0. 


(3.1) 

(3.2) 


This system (3.1)-(3.2) is precisely the NS-Voigt regularization except for the D 1 term included 


in the regularization, and thus we could refer to (3.1)-(3.2) as the NS-Voigt-deconvolution regular¬ 
ization. Since the usual regularization in NS-Voigt is given by just F -1 , and since D « F _1 , this 
“new” regularization D~ 1 F~ l will provide a higher order of consistency for (3.1)-(3.2) to the NSE 


(which are recovered if no regularization is applied). 

In this section we enlarge the derivation of the Voigt models as given in m, and in particular 
we highlight the fact that Voigt models can be seen as ADMs. The Euler-Voigt equations read as 


Vt — a 2 Av t + v ■ Vu + Vq = 
V • Dv = 


(3.3) 

(3.4) 
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and have been introduced in the LES community by Cao, Lunasin, and Titi [13] . They started from 
the inviscid “simplified Bardina”, which is the same model studied in Layton and Lewandowski [55] 
(zeroth order Stolz and Adams) in the case v = 0. The exact words are “... Therefore, we pro¬ 
pose the inviscid simplified Bardina model as regularization of the 3D Euler equations that could 
be implemented in numerical computations of three dimensional inviscid flows. ” They then re¬ 


introduced the viscosity, obtaining the NS-Voigt equations (1.3)-(1.4). Later, in Larios and Titi [341 . 


the model is justified as follows “... it was noted in Cao Lunasin and Titi m that formally setting 
v = 0 in simplified Bardina amounts to simply adding the term —a 2 Au t to the left-hand side of the 
Euler equations yielding what they call the Euler-Voigt equations. Remarkably , if one reintroduces 
the viscous term —vAu the resulting equations happen to coincide with equations governing certain 
visco-elastic fluids known as Kelvin-Voiqt fluids, which were first introduced and studied by A.P. 
Oskolkov J2f 

We show now that the NS-Voigt model can also be derived with usual modeling techniques and in 
the unified framework of deconvolution models. In addition to the motivation of the order reduction 
explained in {25], we follow now a pure approximate deconvolution derivation. 

In order to the derive the model, apply the filter F =: A -1 = (I — a 2 A)” 1 to the incompressible 


NSE (1.11 to get 


A~ 


u t 


+ A -1 [(u • V) u\ - i/A~ 1 Au + VA"V = A" 1 /- 


Rewriting the equation only in terms of the filtered variable (u, p) we get 

u t + A~ x [(Au ■ V) Au] — vA~ 1 AAu + Vp = A -1 /, 

which is the starting point in the Adams and Stolz mm approach. 

Introduce now a family of abstract deconvolution operators {Dm}n which are of order zero 
(explicit examples are those recalled in Section [2A| see especially properties from Lemma 2.1l such 
that 

-Do = I and D m —^ A as A —^ -f-oo ■ 

We then insert the deconvolution into the equations, which allows us to write 

Au ss Dnu. 

Taking the zeroth order approximation, we get (in the new approximate variables v ~ u = D 0 u and 

r~p), 

v t + A -1 [(v • V) v\ — vA~ l Av + Vr = A~ 1 f. 

Finally, by applying the operator A to both sides, we arrive at the NS-Voigt equations 

v t + orAvt + (v ■ V) v — vAv + S7q = /, 


V • v = 0, 


(3.5) 

(3.6) 


Following the same procedure, but using a higher accuracy operator Dm (as in (3.1)) we obtain the 
NS-Voigt deconvolution model 

v t + A -1 [( Dnv ■ V) Dmv\ — oA~ 1 ADmv + X7q = 0. 


After applying A, we recover exactly the RADM model (1.3) studied in Esm]. 


4 Energy analysis 

This section studies the treatment of energy of the RADM system. We begin with the energy balance, 
which is useful both in obtaining a priori estimates for existence theorems and as a starting point to 
study the RADM energy cascade. Then, a Kraichnan-type analysis of the RADM energy cascade is 
given, followed by some numerical experiments to test the energy cascade in forced, homogeneous, 
and isotropic turbulent flows. 
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4.1 Energy balance 


Energy is a quantity of fundamental importance in fluid flow modeling, both in terms of the physical 
relevance of its solutions and its potential for success in simulations. For a model to be successful 
in accurate coarse mesh simulations, it must dissipate sufficient energy to remain stable. We show 
here that the deconvolution in the RADM acts as an energy sponge, which seems to explain a 
smoothing effect from increasing deconvolution in the RADM that was observed in [25] with van 
Cittert approximate deconvolution, and in 1 25] with multiscale deconvolution. 

We present now some integral identities and estimates, which -at the level of Faedo-Galerkin 
approximations- are the core to obtain existence of weak solutions. At the level of weak solutions 
(which are smooth enough) they are still valid and show the balance of the main flow quantities. 


Remark 4.1. Equation (2.11 shows that the RADM conserves a model energy 


E M (t) := ^ (a 2 X7D 1/2 v{t) + D 1/2 v{t) ^ . 


(4.1) 


The conservation of the energy is one of the outstanding open problems for the Navier-Stokes equa¬ 
tions. The analysis of the conservation of the model energy for LES models started with the analysis 
in mu. zy, and several ADM models share this property. 


Remark 4.2. We note that since ||D|| > 1, 


the balance (2.1) suggests that deconvolution in the 

which 


RADM increases the effective viscosity. This is consistent with numerical results in 
reported a smoothing effect by increasing the deconvolution order. 


4.2 Energy spectra 

This section presents an analytic and computational study of the RADM energy spectra. A widely 
accepted theory of turbulent flows starting with intuitions of L.F. Richardson and A.N. Kolmogorov 
is that energy is input at large scales, transferred from large to small scales through the inertial 
range, and finally dissipated at small scales by viscosity [23, 22, 9 . It is further believed that 
viscosity has no effect except on very small scales, on which it acts to decay exponentially fast. 

To fully resolve a turbulent flow, one must resolve the NSE on this entire range of active scales, 
which is far beyond the limit of modern computational tools for high Reynolds number problems. 
Computing the NSE without full resolution is widely known to be unstable. For a model to be 
computable in practice, its range of active scales needs to be significantly smaller than for the NSE, 
and thus studying its energy spectra can help make such a determination. 

Since the RADM takes a form similar (in some sense) to that of Leray-a, NS-a, and the simplified 
Bardina models, we adapt physics-based energy spectra analysis performed in ns in hd for these 
respective models. We begin with a decomposition of the velocity into its Fourier modes, which 
gives the energy balance 

\^t {( Vk ’D N (k) v k) + a 2 {-Av kl D N {k)v k )^J + v(-AD N (k)v k ,D N (k)v k ) =T k -T 2k , (4.2) 

where D]\r(k) is the Fourier transform of the deconvolution operator Dn, and 


Tk ■= -((D N (k)v£ ■ S7D N (k)v k , D N (k)v k ) + {(D N (k)(v k + vjf ) • S7D N (k)(v k + vjf ), D N (k)v ^), 


with 

v k '■= V T and v k '■= J2 V T 

\j\<k \j\>2k 

Note that T k — T 2k represents the net amount of energy transferred into wave numbers between 
[fc, 2k). Time averaging the energy balance (4.2) gives the energy transfer equation 


(i/(-A D N (k)v k ,D N (k)v k )} = (T k ) - (T 2k ). 


(4.3) 
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Since the RADM model energy is defined by (4.1), as this is the conserved quantity in the absence 
of viscosity and external forces, we take D = D jy and define the model energy of an eddy of size 
k~ x by 

E M (k) ~D N (k)(l + a 2 k 2 ) J2 N 2 - 

lil=fc 


Combining this with the energy transfer equation (4.3) gives the relation 


(T fc ) - (r 2fc ) 


vk 3 Djs[[k)EM{k) vk 3 (N + 1 ).Em(&) 


1 + a 2 k 2 


1 + a 2 k 2 


with the last relation coming from Lemma 2.1 since D^ik) ~ (N + 1) for sufficiently large k. If the 
wave number k belongs to the inertial range, then (Tk) ~ (TA) and there is no leakage of energy 
through dissipation. Then, we have that k is in the inertial range if 

uk 3 (N + l)E M (k) 


1 + a 2 k 2 


< (T fc ). 


This shows that increasing the deconvolution order N has the effect of shortening the inertial range 
of the RADM. This is consistent with the energy balance given above, which suggests that increasing 
N reduces the total energy and increases the effective viscosity. 

To determine the kinetic energy distribution in the inertial range, we begin by defining the 
average velocity of a size fc -1 eddy, and relating it to model energy via 


U k ■= (( v k ,v k )) 1/2 


’ f 2k E M {k) \ V2 

Jk (1 + a 2 fc 2 )Zbv(fc) J 


f kE M (k ) \ 1/2 

\ (1 + a 2 k 2 )(N + 1) / 


The total energy dissipation rate is determined by the energy balance to be 


e M := WIVTMI 2 ), 


and from the Kraichnan energy cascade theory gives that the corresponding turnover time for eddies 
of this size is 

1 _ (l + fc 2 a 2 ) 1 / 2 (A^+l) 1 / 2 
Tk ''~kU~ k - k 3 / 2 E M {ky/ 2 ’ 

the model energy dissipation rate Em scales like 


1 r 2k 

£m ~ — / E M {k) 
T k Jk 


k 5 / 2 E M (k) 3 / 2 
(l + k 2 a 2 ) 3 / 2 {N +1) 1 / 2 ' 


Solving for EM(k) provides the relation 


2/3 

e M 


(1 


E M (k) 

which implies a kinetic energy (E = 4 ||f|| z ) spectrum 


a 2 k 2 y/ 3 (N ■ 
fcs/ 3 


l) 1 /3 


E(k) 


E M {k) 

(1 + a 2 k 2 )D N (k) 


2/3 

e M 


fc 5 / 3 (l + a 2 k 2 ) 2 / 3 (N + l) 2 / 3 


As is the case for related models, such as those of a-type or the simplified Bardina model, the kinetic 
energy spectrum can be divided into 2 parts. If ka > 0(1), then 


m 


2/3 

e M 


2/3 

e M 


fc 5/3( Q ,2 fc 2)2/3( iV + !)2/3 k 3 CX 2 (N + l) 2 / 3 ’ 

















and if ka < 0(1), then we get 


2/3 

E(k) ~ _ — _ 

1 ’ fc 5 / 3 (TV + l) 2 / 3 ' 

These scalings suggest that for larger scales in the inertial range, a k ~ 5//3 roll-off of energy is expected, 
but for smaller wave numbers, the slope increases to k~ 3 . This implies that significantly less energy 
will be held in higher wave numbers, suggesting the model is indeed more easily computable than 
the NSE, which has a fc -5 / 3 roll-off of energy through its entire inertial range. 

4.3 Numerical testing of energy spectra 

To test some of these scaling results, a direct numerical simulation (DNS) of the non-dimensionalized 
RADM is performed for forced turbulence, which is one of the most extensively simulated turbulent 
flows (see e.g., d Hi 1301132 HZj). In particular, we study the energy spectrum of the model for 
three-dimensional homogeneous and isotropic turbulent flow, on a periodic cubic box with a side 
length 27r, with Reynolds number Re = 200. 

We use a pseudo-spectral method for the spatial discretization of the model with full de-aliasing, 
and a second-order Adams-Bashforth scheme for the time stepping. In addition, we employ a forcing 
method used by Chen et al. mm and She et al. [43 j. The numerical forcing of a turbulent flow is 
the artificial addition of energy at the large scales in a numerical simulation. Statistical equilibrium 
is achieved by the balance between the input of kinetic energy through the forcing and its output 
through the viscous dissipation. Van Cittert deconvolution is used with order N = 0,1,2, and we 
will specify TV in the plots. 

We first begin by testing the scaling law in the previous section, which predicts that the kinetic 
energy scales with fc -5 / 3 at the beginning of the inertial range, and transitions to k ~ 3 near the end 
of the inertial range. A plot of the resulting energy spectra from a computation with TV = 2 and 
a = 1/32 at the resolution of 128 3 is shown in Figure[l] Along with the spectrum, on the log-log plot 
are also lines with slopes -5/3 and -3, and we observe that the spectrum appears to be in agreement 
with these slopes near the beginning and end of the inertial range, respectively. 



Figure 1: The energy spectrum of the RADM at 128 3 resolution (TV = 2, a = 1/32). Slopes on the 
log-log plot appear in agreement with the predicted slopes of —5/3 and —3 at the beginning and 
end of the inertial range, respectively. 


We also investigate the sensitivity of TV on the energy spectra of the RADM. Here, we compute 
with 128 3 resolution, a = 1/16 and TV = 0,1, 2. Results are shown in Figure [2j and we observe (as 
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our analysis predicts) the RADM energy curves shorten as N increases; that is, as N is increased, 
there is less (or equal) energy in each wave number across the entire spectrum. This shows there 
is less total energy in the N = 2 system than N = 1, and N = 1 compared to N = 0, which is in 
agreement with the above analysis. 



Figure 2: Energy spectra for values of N = 0, 1, 2 at 128 3 resolution for a = 1/16, along with 
the energy spectrum obtained from the DNS of the NSE at 256 3 resolution. The proposed model’s 
energy spectra become closer to the DNS energy spectrum with increasing N. 

Lastly, we use this experiment setup to test the accuracy of the RADM energy spectra against a 
DNS of the NSE. In Figure[3j we show the resulting energy spectra for the RADM with various values 
of a for N = 2 at the resolutions of 63 3 and 128 3 along with the DNS energy spectrum obtained 
from the NSE at 256 3 resolution. The energy spectra are provided for three cases of a = 1/8, 1/16, 
and 1/32. As a decreases, the energy spectra of the RADM become closer to the highly resolved 
DNS energy spectrum. 


5 Existence of global attractors for the RADM 


The long-time behavior of evolution equations can be characterized by their attractors, which in 
many important cases are finite dimensional. In particular, the notion of global attractor (cf. [46] ) 
seems the most relevant in the context of the NSE. It is well known that such an attractor exists for 
the 2D NSE, while the lack of well-posedness in the 3D case is reflected in partial results concerning 
weak or trajectory attractors. In the case of LES models, the solution has naturally better regularity 
properties (even if a certain hyperbolic behavior is expected), hence it is natural to expect existence 
of such an object. In particular, for the NS-Voigt model the analysis of these questions has been 
addressed in [22J [25] , where existence and regularity of the global attractor and also determining 
modes are studied. Since the NS-Voigt model is a peculiar system (pseudo-parabolic, see also 


results in Section 6.2) with also non-viscous features, the proof of existence of the global attractor 


uses techniques typical of hyperbolic systems and the semigroup results to be only asymptotically 
compact. In particular see the questions related with the regularity of the solution in the space 
variables recalled in Sec 16.11 

In this section we follow very closely the approach in (2® and we just highlight the changes 
needed to adapt the proof to the new RADM system. The main difference with respect to the latter 
reference is that, due to the presence of approximate deconvolution operators, we need to deal with 
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Figure 3: Energy spectra for various values of a = 1/8, 1/16, and 1/32 for N = 2 at (a) 64 3 and 
(b) 128 3 resolutions along with the energy spectrum obtained from the DNS of the NSE at 256 3 
resolution. The proposed model’s energy spectra become closer to the DNS energy spectrum with 
decreasing a. 


the space-periodic case. 


For simplicity of analysis and notation, following [2S] in this section we consider the model (1.3)- 

(5.1) 


(1.4) in the equivalent form 

v t + a 2 Av t + P(Dv ■ X7Dv) + vADv = Pf, 


where P is the Helmholtz-Leray orthogonal projection in L 2 (T) onto Vo, and A := —P A is the 
Stokes operator with periodic boundary conditions. The existence of a global attractor will follow 
from a well-known result, for which we refer for instance to H5i . 

Theorem 5.1. Assume a semigroup S(t) : V 1 V 1 , 0 < to < t can be decomposed as 


S(t)=Y(t) + Z(t), 

where Z(t) is compact in V 1 for each 0 < to < t. Assume also there exists a continuous k : 
[to, oo) x R + —> R + such that for every r > 0, k{t, r) —> 0 as t —> oo and that 

\\Y(t)v\\ v < k{t,r), Vt > t 0 , Vw € V 1 : ||u|| < r. 

Then S(t) : V 1 —> V 1 , t > 0, is asymptotically compact. 


The main result of this section is the following. 


Theorem 5.2. 

Let S(t) : V 1 


Let vq 
V 1 , 


G V 1 

t G 


and assume the forcing is only spatially dependent, Pf = Pf(x ) G V . 
R + be the continuous semigroup generated by the problem (5.1) (the 


existence of which is guaranteed by Theorem 2.1 of m- Then S(t) has a global attractor. 


Proof. As usual we begin by showing the existence of an absorbing ball. Then, we show that the 
assumptions of Theorem |5.1| are satisfied, to conclude that a global attractor exists. We proceed 
formally, but calculations again can be completely justified by a Faedo-Galerkin approximation. 
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We multiply (5.1) by Dv and proceed as above for the energy inequality analysis, along with 
Cauchy-Schwarz, and Young inequalities on the forcing term to obtain the bound 


dt 1 01 


\7D 1/2 v 




+ v\\S7Dv\\ 2 < -\\Pft,. 


Using that D is positive and self-adjoint, along with assumptions on its boundedness, gives 


dt “ 


\7D 1/2 i 




v do 


VD 1/2 ? 


2 1 9 

<-|ra_r, 

v 


and then applying Poincare provides the bound 


dt ' “ 


Vd 1/2 i 


D x ' 2 x 


2 \ / a 2 


Setting co = min{a 2 , A}, we get 


d 

dt 


\7D 1/2 v + D 1/2 v 

and thus from the Gronwall inequality, 


+ vdn 


+ Cq [ a 


Vfl 1/2 i 


+ A 


D l ' 2 x 


\7D 1/2 v + D 1/2 i 


<r\\ p f\\-i- 


< -\\PfW-! 

V 


limsup fa 2 VD 1/2 v(t) \ D 1/2 v(t ) ) 

£-H-oo \ ) 


< 


\\PfW-! 

uc 0 


Again applying the lower boundedness assumption on D, this reduces to 

2 \ 2nrax{a 2 ,A -1 } 


limsup (a 2 ||Vv(f)|| 2 + ||v(f)||' 

£—>■00 ' 


< 


» 2 dl 


Pf ll-i 


which implies that S(t) : V 1 —> V 1 , t G R + has as absorbing ball Bi := 5(0, Mi) C V 1 , and we 
have the uniform estimate 

IK*)lli < M u 

for t large enough. 

We now show the assumptions of Theorem 5.1 are satisfied, following analysis similar to [251 . 
Observe that S(t) can be represented as 

S(t)v o = Y(t)v 0 + Z(t)v 0 , 

where Y ( t ) is the semigroup generated by the linear (which is linear since also D is a linear operator 
when the van Cittert deconvolution is used) problem 


y t + a Ay t + vADy = 0, 

2/(0) = vo, 


while z(t ) = Z(t)v o solves 


z t + a 2 Az t + vADz = Pf — P(Dv ■ VDv), 

3(0) = 0, 


(5.2) 

(5.3) 


(5.4) 

(5.5) 


with v solution of (|5.1|). We take the L 2 -inner product of (|5.2|) with y to find 

2 


d 

dt 


2 + a 2 ||Vy(f)|| 2 ) + 2v |V5 1//2 y 


= 0, 


and then using that D is positive and bounded from below, along with the Poincare inequality, we 
obtain as in the previous calculations, 


d 

dt 


(||y(t)|| 2 + a 2 |! Vy(t)|| 2 ) + c 0 (A ||y(f)|| 2 + a 2 ||Vy(t)|| 2 ) < 0, 


12 








































































which immediately implies 


(||j/(t)|| 2 + a 2 I!Vy(t)|| 2 ) < e~ Cot (|K|| 2 + a 2 ||Vz, 0 || 2 ) . 


Therefore the semigroup Y ( t ) is exponentially contractive. 

By using the properties of D, and standard inequalities in Sobolev spaces (cf. [551 Eq. (3.11)]) it 
follows that 


\\P(Dv ■ VUu)||y-l/2 


<t> ev 1 


sup 

A\au*4,\\=i 


\v • VDu, P' 


<C\\Dv\\l4A^\\ 

< Cdl\\v\\yl || A l/4 (/)||, 


hence showing that, since v is a weak solution to (1.31, then 


P{Dv ■ VDv) G L°°(0,T; V 1/2 ), 


for all positive T, and with bounds independent of T, but depending essentially on ||D||. 

Then we take (5.41-( 5.5) and we test it by A 3 / 2 z (procedure which can be again justified by 
approximation with Faedo-Galerkin method). We observe that Mv = (A‘z,z) = \\A'/*zf, by 
definition. Hence, by noting that Pf — P(Dv ■ VDv) G L°°(0, T; V 1 / 2 ) we get 


\j t (o 2 ||A 3/1 2|| 2 + M 1/4 *l| 2 ) + v\\A 3,i D l,2 z\\ 1 

= - ^ (« 2 ||s|lv.3/2 + llsllv-i/s) + I'll-D 1,2 z|l?,.„ 

= (Pf - P(Dv ■ VDv), A 1/2 z) 

= ( A~ 1/4 {Pf - P(Dv ■ VDv)),A 3/4 z) 

< (\\f\\v-v* + II P(Dv ■ VI>t;))|| v - 1/3 )[|*||v3/ a . 


Using previous estimates we obtain 


^ ( a “ll 2 llv3/2 


Wv 1 / 2 ) 


/HU 1 / 2 


Z U v 3 /- 




1/2 


dj | 


llv 1 ) 


and the bound on Iklk 1 < Mi implies that 2 is bounded in U 3 / 2 . This means that Z(t) maps V 1 into 
U 3 / 2 , and thanks to the compact embedding U 3 / 2 c U 1 , Z(t) is compact. Thus the assumptions 
of Theorem 5.1 are satisfied, and we conclude that S(t) is asymptotically compact. By J35], we 
conclude that S(t) has a global attractor in V 1 ; the rest of the proof follows as in [25p With the 
same arguments one can also show that if the force is smoother, say in Pf G V 1 , then the attractor 
is a compact set in U 2 . □ 


6 Pipe flows 

In this section we perform some analytical computations with very classical tools, in order to focus 
on some special analytical features of the solutions of Voigt models and also to better understand the 
possible role of the RADM (and thus the NS-Voigt system for Dg = I ) in modeling of visco-elastic 
flows. We recall that prior to turbulence modeling, NS-Voigt equations appeared as a linear model 
for certain non-Newtonian fluids, see a series of papers by Oskolkov [3B] [32] 

In order to understand, at least qualitatively, the effect of the (non-viscous) damping term 
—a 2 A u t on the behavior of solutions, we restrict to the very simple setting of: a) fully developed 
flows, that is the velocity is directed only along the axis of the pipe; and b) we restrict to the NS- 
Voigt case; that is, assume D = D 0 = I. The deconvolution order N > 0 could be used, but details 
would change for each specific deconvolution operator, and also the analysis will become (even more) 
technical and also boundary conditions for filtering at the wall would need to be specified. Under 
the two above assumptions we develop analytical solutions, and even if we are aware that these 
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solutions are not turbulent and describe only the linear behavior, they constitute a good source 
of data for debugging complex codes, as only very few analytical solutions are known to exist. 
Moreover, analytical solutions can be used to gain some physical insight into the model’s solutions. 
On the other hand, we recall that a remarkable feature of the NS-Voigt equations is that is one 
of the few LES models which can be fully treated with Dirichlet boundary conditions and with 
given external force, see nans]. Since the NS-Voigt equations can be treated in the presence of 
solid boundaries with Dirichlet boundary conditions, it is natural to investigate some class of exact 
solutions for NS-Voigt, as in the case of the NSE, to directly compare the results. 

6.1 Remarks on the regularity of the solution 

In this section we make some remarks on the solution of the NS Voigt equations, but similar reasoning 
can be transferred to the RADM (at least in the space periodic-setting, while here we have a problem 
with solid boundaries and zero Dirichlet conditions). In particular, as remarked in Larios [33j the 
NS-Voigt have a pseudo-parabolic behavior, as the abstract equations considered in Carroll and 
Showalter m- One of the most interesting point is that the solution of the NS-Voigt preserve the 
same regularity as the initial datum. Contrary to the NSE for the solution of the NS-Voigt it is not 
possible to show an improvement of the regularity, especially for what concerns the space-variables. 
We recall that for the NSE (or more generally parabolic equations), as soon as a weak solution is 
strong, then it becomes smooth (even C°° if the force and the boundary of the domain are so) for 
all positive times. This has particular consequences on the fine properties of weak solution of the 
NSE, when obtained as limits of solutions to the NS-Voigt equations as a —> 0 + , see JTTj . 

Here we show with a particular simple example that the nonlinear NS-Voigt cannot produce an 
improvement of the regularity, by making some very explicit computations on the spectrum of the 
Laplace operator. This can be also used to shed some light in the different behavior of the NS-Voigt 
with respect to the NSE, concerning some very special properties. 

Let us consider the NS-Voigt system in a special physical situation. We have a channel with 
cross section E and with axis directed along the * 3 -direction. We consider an incompressible fluid, 
modeled by NS-Voigt equations in a semi-infinite straight pipe P = E x R + cl 3 . In a reference 
frame with z (we call z := £3 for historical reasons) directed along the pipe axis and x := (£ 1 , 2 : 2 ) 
belonging to an orthogonal plane, the NS-Voigt equations read 

v t — a 2 Av t + (v ■ V) v — uAv + Vp = / (x, z) £ E x M + , t £ R + , 

V ■ v = 0 (x,z) £ E x K+, t £ K+, (6.1) 

v = 0 (x, z) £ dE x R + , t £ R + , 

where v(t, x, z) and p(t, x, z) respectively denote velocity and pressure, and v > 0 represents kine¬ 
matic viscosity. 

Remark 6.1. In this subsection and later on also in Section 
remember that the velocity is a vector field and that we will look 
velocity w(t,x,z) is the only non-zero one. 

By following classical calculations dating back to Hagen and Poiseuille, we look for fully-developed 
solutions (also named Poiseuille-type solutions): 

v(t, x, z) = (0,0, w(t, x)) and p(t, x, z) = — A(i, x, z) + po(t), 

where po(t) denotes an arbitrary function of time. The Poiseuille-type ansatz implies that the 
convective term cancels out identically, and this motivates the statement that the results will concern 
only the linear behavior of solutions. Moreover, it is easy to deduce from the equations that the 
pressure is p(i, z) = —A (t.) z. Finally, the dependence of w on the space variables X\ and £ 2 allows 
us to consider a problem reduced to the cross-section E : In the so-called “direct problem,” without 
external force, a given pressure-drop is assigned and the problem is the following: Given A (t), find 


6.2 we use the symbol v(t,x,z) to 


for special solutions, where the axial 
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w(t, x) such that 


x £ E, t £ K+, 
x e &E, f e R+ 
x £ E, 


( 6 . 2 ) 


! w t (t, x) - a 2 A x w t (t, x) - uA x w(t, x) = \{t), 
w(t, x) = 0 
u;(0,a;) = w 0 (x) 


where A x denotes the Laplacian with respect only to the variables x\ and cc 2 - 

In this particular setting the nonlinear term is identically vanishing, but we still have a boundary 
value problem. We will start by making some explicit calculations in the case A (t) is extremely 
smooth, say A = 1 and we assign an initial datum wq(x) £ Hq(E). This corresponds to the classical 
pressure drop problem with a constant gradient superimposed to sustain the flow, but as we will see 
this assumption is completely inessential. We have the following theorem. 


Theorem 6.1. Let E CM 2 be smooth and bounded. Let v(0,x,z) = (0,0, wq(x)) be given such that 
Wq(x) £ Hq(E), and let A = 1. Then, there exists a unique solution to (6.2) such that 


w £ L°°(Q,T-Hl{E)) VT> 0 


hence a unique solution to the problem (6.1) under the Hagen-Poiseuille ansatz. 
u>o ^ H S (E) for any s > 1, then also w[t,x ) ^ H S (E) for all t > 0. 


Furthermore, if 


The lack of regularization in the NS-Voigt equation is due to the special role of the inviscid 
regularization and it has some consequences when the Voigt model is used for theoretical purposes, 
see for instance m- 

of Thm. \ 6.1\ First we note that the existence and uniqueness of solution with initial datum wq and 
external force / = (0,0,0) follows as in [13] . Here we make a more explicit construction in order 
to highlight the regularity properties. The proof is based on an explicit analysis of the frequency 
budget, working with a spectral basis to construct the solution. (In the space-periodic case this can 
be done even more explicit with a Fourier series expansion, but here we deal with the boundary 
value problem). 

Let {< f > m } denote the eigenfunctions of the Laplace operator with Dirichlet conditions on E, 


A(fm — A rnf^m hd, 

(j) m = 0 on dE , 


and by standard spectral theory the scalar functions {<t> m } ar e smooth, orthonormal with respect to 
the L 2 (E) scalar product, and the eigenvalues {A m } are an increasing sequence of strictly positive 
terms. Hence we construct the solution as w(t, x) = J2 m >i c m{t)(j)m{x) and it is immediate to check 
that we have to solve the following ordinary differential equation, for all m > 1 

(1 + a 2 A m )c , m (t) + z'c m (f) = /3 m , (6.3) 

where fd m = f E 4> m {x) dx and with the initial condition c m (0) = f E Wo(x)<fi m (x)dx. Observe that, 
from the regularity hypotheses on the initial datum, we have 

y: A m |c m (0)| 2 < +oo A^|c m (0)| 2 = +oo, for all s > 1. 

ra>1 m> 1 


Explicit integration of (6.3) gives 


c m (f)=c m (0)e !+“ 2A m + j3 n 


m > 1. 


It is straightforward to check that X)i< m <A r c m(t)(f> m (x) —> w(t,x) as N —► +oo, and from the 
explicit expression of the solution we can see that even if /?„, would be zero, then ||u>(f, x)||jjs = oo, 
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_ A m v -f- 

for all s > 1, and for all t > 0, due to the fact that the exponential e >+" 2> ™ does not produce any 
smoothing (as if a = 0), since 

Mm < A m V < V_ 

1 + a 2 Xi ~ 1 + a 2 X m ~ a 2 

The situation will be different for what concerns time regularity, since any time-differentiation will 
produce the zeroth order term w hich will not affect the convergence of the series describing 

the solution. □ 


6.2 Pulsatile flows 


In this section we consider pulsatile flows, which are driven by a time-periodic force (the pressure 
gradient). This is a remarkable case motivated also by recent studies in 0||B] for heart-beat driven, 
human physiological flows, including blood and cerebrospinal fluid. 


In the so-called ‘direct problem,’ a pulsatile pressure-drop is assigned in problem (6.2). Together 


with this problem, in many cases of physiological interest (when the measure of the pressure cannot 
be done directly with sufficient accuracy by phase-contrast Magnetic Resonance Imaging or Doppler 
ultrasonography) it will be also meaningful to study a problem with assigned flux: f E w(t, x) dx = 
<3>(£), for some smooth-enough given scalar function Q(t). The Poiseuille-type ansatz will lead to the 
problem: Given $(f), find (w(t,x),X(t)) such that 

'd t w(t, x) — uA x w(t,x) = X(t), x £ E, t £ R + , 

w(t,x) = 0, x £ dE, t £ R + , 

^f E w(t,x)dx = &(t), i£l + , 

which is an inverse problem and it is linked to one of the nowadays classical Leray problems. 

The study of this problem when E is a circle, A (t) is time-periodic, and a = 0 dates back to 
Richardson (1926) and Sexl (1930), while a big interest for applications to bio-fluids raised especially 
after the work of the physiologist J. R. Womersley S3- In particular, he quantitatively clarified why 
for the NSE, if the pulsation is fast enough, then the qualitative behavior of solution can change 
drastically: Instead of a parabolic type profile (as in the stationary case) a new phenomenon called 
“Annular effect ” appears: If the frequency of pulsation of A is large enough (with respect to the 
other relevant parameters of the problem), then the maximum of the velocity is not located along 
the axis of the channel, but near the wall. 

This is particularly of interest in blood modeling, where -at first approximation- we have a 
pulsatile flow (blood driven by heart beat) in a straight channel (larger arteries). This motivated 
the analysis of Womersley, who explained by analytical formulas the experimental observation that 
in larger arteries of the rabbit and the dog there is also a reversal of the flow. This suggests that 
pulsatile flows, even in the linear range, can be different and much more complicated than the 
steady ones. The analysis of Womersley was based on a single mode pressure-drop A(t) = e mt and a 
circular cross-section E = {x\ + < R}. In cylindrical coordinates (r, 9) and after the separation 


of variables w(t,x) = e lnt W(r), the system (6.2) becomes: find W{r) such that 


inW(r) — v(w”(r) 


+ 


W{r) 


= 1 


for 0 < r < R, 


W\ 0) = W{R) = 0. 


The above equation for a cylindrical symmetric solution can be solved by means of an expression 
involving the zeroth order Bessel function Jo(r). The most widely known contributiorj^jof Womersley 
is that of determining the non-dimensional quantity (named later Womersley number) 


W„:=R^, 

^^Most important to this field, but recall that he recruited A.M. Turing for the code-breaking project in WWII. 
He also convinced Sir G. Darwin to start the project of an all British computer and coined the name Automatic 
Computing Engine (ACE) for an early electronic computer. 
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which determines the qualitative behavior of solutions, where R is the radius of the channel, ui the 
frequency of the pressure drop, and v the kinematic viscosity. 

The Womersley number is the ratio of the transient or oscillatory inertia force to the shear force. 
When Wo < 1, the frequency of pulsations is sufficiently low that a parabolic velocity profile has time 
to develop during each cycle, and the flow will be very nearly in phase with the pressure gradient. In 
this case using instantaneous pressure to modulate a parabolic profile will give a reasonable solution. 
When Wo > 10 the frequency of pulsations is sufficiently large that the velocity profile is relatively 
flat or plug-like, and the mean flow lags the pressure gradient by about 90 degrees. 

In addition to this peculiar effect, the Womersley solution has become a paradigm in the analysis 
of biological fluids (see Fung [2J and Quarteroni and Formaggia uni) for its simplicity, being one 
of the very few exact unsteady solutions to the NSE, to debug complicated CFD flow, and also to 
provide an improved source of initial/boundary data, as well as a benchmark solution for pulsatile 
flows. 

We now perform an analysis of the NS-Voigt equations for this problem. Assume that = e lut , 
and by the ansatz w(t,x ) = e lut W{r) in the case of (1.3) we arrive at the following boundary value 
problem for a single ODE: Find W(r) such that 


W (r) — a 2 (w"(r) 


W(r) 


'(W"(r) 


W (r) 


= 1 


for 0 < r < R, 


W'( 0) = W(R) = 0. 


The above equation can be exactly solved by means of an expression involving the zeroth order 
Bessel function J 0 (r). By explicit calculations, we get 

w(r) = — i— 

lw \ 

\ \ VC* M-llS , 


which reduces to the classical solution obtained by Womersley in the case a = 0. Hence, we get that 
the natural non-dimensional quantity is the following (which we call a-Womersley number) 


a-Wo := R 


\J a 4 u > 2 + v 2 


The a-Womersley number is always smaller than the Womersley number of the flow with the same 
viscosity, radius of the pipe, and frequency. In fact 


Wo iT a 4 uj 2 

a-Wo V v 2 

and the difference can be significant if a is not properly chosen. This gives the hint that a not finely 
tuned choice of a can completely change the behavior of solutions, since the flow can pass from one 
regime to another, if Wo is large and a is not small enough. 

Here, we want to give a more explicit and real solution by considering the 2D case, that is, the 
cross section E is the ID interval E =:] — R,R [, and showing the effect of a pulsatile force cos (uit), 
which is not a peculiar effect of the space dimension. In particular, we will see that the qualitative 
behavior is not better in a 2D pipe, with respect to the 3D one. We consider then the following toy 
problem 


I d t w(t, x) — a 2 w tX x(t, x) — vw xx {t, x) = cos(wf), ire]— R, i?[, t € R + , 

| w(t,x) =0 x = ±R , t € K + , 

which corresponds to the pulsatile flow in a rectangular 2D channel. In this case the exact solution 
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(which is real as can be checked by direct inspection) has the following explicit expression: 


. . sin (u() 1 

w{t, x) = -1- — (—sm(u;f) — i cos(w t)) sech 

uj 2c u 


R 


cosh 


d-(— sin(fu) + i cos(w t)) sech 

2uj 


R 


cosh 


If we set, in the case a = 0 (i.e. the NSE), u = 1, R = 1, and w = 144, we find Wo=12, and the 
solution at t = 0 shows the annular effect and reversal of the flow at distance about R/ 2, as shown 
in Figure [4j 



Figure 4: Shown above is the profile of the Navier-Stokes solution at time t = 0 with Wo = 12. 

Next, we plot with the same parameters, the corresponding solutions in terms of different values 
of a. Pictures in Figure [672] show that the qualitative behavior is quite different with respect to the 






Figure 5: Profile of the solution at time t = 0 with (a) a = 1, hence a-Wo = 0.999988; (b) a = 1, 
hence a-Wo = 1.99961; (c) a = 0.1, hence a-Wo = 9.06295 and still no reversal starts (even if the 
solution at x = 0 is about 10~ 6 ); (d) a = 0.1, hence a-Wo = 11.6399 

variations of a. Simple calculations show that, in order to reconstruct the correct pulsatile behavior 
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with flow reversal, the parameter a must satisfy 


y^Wo 4 - 10000 

a < 

As is evidenced in Figure |th2| this puts strong limitations on the choice of a, and moreover such a 
small a moves the NS-Voigt systems closer and closer to the NSE. Moreover, in LES, a has the role 
of the smallest resolved length scale, and when looking at the system in this sense, the limitations for 
this model are those of a good resolution of the flow, otherwise transient features are lost. This could 
be used to speculate that the space filter alone cannot be very good when studying time-evolution 
problems with transient or periodic effects, or at the very least the smallest resolved length scale 
should not be constant across the domain. This idea is mentioned in [26], and indeed better results 
for turbulent channel flow were found when a was chosen of the order of the local mesh width, which 
was much finer in the boundary layer. 


7 Conclusions 

We have presented a detailed mathematical analysis of several important aspects of the RADM. 
In particular, we presented a detailed analysis of the model’s treatment of energy, which showed 
that deconvolution acts to drain energy from the system, and likewise affects the energy spectrum. 
We derived a scaling for energy in its Fourier modes, which shows that increasing deconvolution 
in the model shortens the inertial range, and that the inertial range is split into two parts: the 
larger scales of the inertial range have a fc -5 / 3 scaling with energy, while the lower numbers have a 
faster & -3 scaling. Interestingly, this is the same spectral scaling as both the Leray-a and simplified 
Bardina models. We presented energy spectra from computations of forced, isotropic, homogeneous 
turbulence with varying N and a, which verified the predicted scalings, and showed good agreement 
on large scales with DNS of the NSE. In addition to the RADM energy study, we also proved 
existence of smooth global attractors, and provided an analytical study of pulsatile flows that shows 
small a may be requested to study certain flows. 

There are several important directions for future work with the RADM. First and foremost, 
more benchmark testing needs done. In particular, since the model performed well with Re T = 180 
turbulent channel flow, future testing on the benchmarks of Re T =395 and 590 are natural next 
steps. Turbulent flow around obstacles is another important test problem, as is applying the model 
to specific application problems such as flow through an aorta. The model could also be studied 
with different types of filtering and deconvolution. 
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